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Abstract 

We apply the holonomic gradient method to compute the distribution function of a 
weighted sum of independent noncentral chi-square random variables. It is the distribu¬ 
tion function of the squared length of a multivariate normal random vector. We treat this 
distribution as an integral of the normalizing constant of the Fisher-Bingham distribution 
on the unit sphere and make use of the partial differential equations for the Fisher-Bingham 
distribution. 

Keywords and phrases: algebraic statistics, cumulative chi-square distribution, Fisher-Bingham 
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1 Introduction 

The weighted sum of independent chi-square variables appears in many important problems in 
statistics. In the problems for testing against ordered alternatives, cumulative chi-square statistic 
(cf. [7], [13] ) has a good power. For studying the power function of the cumulative chi-square 
statistic, we need to evaluate the distribution function of a sum of weighted independent noncen¬ 
tral chi-square variables. Goodness of fit test statistics based on empirical cumulative distribution 
function, such as the Cramer-von Mises statistic or the Anderson-Darling statistic OH), are infi¬ 
nite sums of weighted independent chi-square variables. Chapter 4 of [4] gives a survey of these 
statistics. Under an alternative hypothesis the chi-square variables are noncentral. For studying 
the power function of these statistics we want to approximate the infinite sum by a finite sum of 
sufficiently many terms and compute the cumulative distribution of the finite sum. 

An exact evaluation of the cumulative distribution function of a weighted sum of independent 
noncentral chi-square random variables was considered to be a difficult numerical problem (see 
0)- Although the moment generating function is explicitly given, its Fourier inversion to evaluate 
the density function and the cumulative distribution function is difficult as extensively discussed 
in Chapter 6 of [IB]. See [3] for the similar problems in other areas of applied mathematics. 

Recently in p3] we proposed the holonomic gradient method (HGM) for calculating distri¬ 
bution functions and the maximum likelihood estimates using differential equations satisfied by 
a probability density function with respect to the parameters. Since then the method has been 
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successfully used in many problems, including the computations related to the Fisher-Bingham 
distribution on the unit sphere ([lUj, [H], [9], [[15j). In this paper we utilize the results on HGM 
for the Fisher-Bingham distribution to evaluate the distribution function of a weighted sum of 
noncentral chi-square random variables. 

Let X denote a d- dimensional random vector following the multivariate normal distribution 
Consider the cumulative distribution function G(r) of ||X||: 

G(r)= H (x ■ '' )Te_1(x ■ "0 dx - (1) 

We call G(r) the ball probability with radius r. By rotation we can assume that E = diag(oy,..., af) 
is a diagonal matrix without loss of generality. Hence G(r) is the distribution function of the square 
root of a weighted sum of independent noncentral chi-square random variables, where weights are 
erf, i = 1 ,d. Furthermore the conditional distribution of X given its length r = ||X|| is the 
Fisher-Bingham distribution. This fact allows us to directly apply the results for the Fisher- 
Bingham distribution to the evaluation of the distribution of the weighted sum of independent 
noncentral chi-square random variables. As we show in Section H] our method works very well, 
both in accuracy and speed. 

The organization of this paper is as follows. In Section [2] we summarize known results on HGM 
for the Fisher-Bingham distribution and show how they can be used to evaluate the distribution 
of the a weighted sum of independent noncentral chi-square random variables. We also discuss 
the problem of initial values needed to use HGM. In Section [3] we present asymptotic results for 
the Fisher-Bingham integral and its derivatives for the case that the length of the multivariate 
normal vector diverges to infinity. This result is used to check the the numerical accuracy of our 
experiments in Section [4l We end the paper with some discussions in Section [5} 

Acknowledgment. This work is supported by JSPS Grant-in-Aid for Scientific Research No. 
25220001 and Grant-in-Aid for JSPS Fellows No. 02603125. 


2 Holonomic system and initial values 

Let 

E = diagOi,..., a d ), h = (AW • • •, dd) T ■ 


We define new parameters A*, t*, i — 1,..., d, by 

X - 1 _ Ah 

2 °f 

and the Fisher-Bingham integral /(A,r, r) by 


/(A,r,r) 


/ d d \ 

/ exp ( V A $ + V TiU ] dt, 

J S d ~Hr) Vi=l i=l / 


( 2 ) 


where A = (Ai,..., Ad), r = (ri,..., r d ), S d x (r) = {t e | t\ +-f t 2 d = r 2 } is the sphere of 

radius r and dt is the volume element of S' d_1 (r) so that 


dt = r d ~ 1 S d -i, S d -! = VoK^-^l)) 


•d-li 


' S d hr) 


27 T d / 2 

w 
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Then G(r) in (JTJ) is written as 


G(r) = 


nt.v^Ai 


7 T 


d/2 


exp - 



T" 1 / f(X,r,s)ds. 


(3) 


We will numerically integrate the right-hand side of (|3]). We denote the partial differential operator 
with respect to A by d\. For t e S' d ^ 1 (r), (t\ + • • • + t 2 d )/r 2 = 1 and 


Titij dt 

= ^(0Ai + -" + dx d )f{X,T,r). (4) 

By HGM we evaluate d\ i f(X, r, r), i = 1,..., d, and use (J3| to compute /(A, r, r). In fact we also 
evaluate d n f( A, r, r), i — 1 ,..., d. 

Define a 2d-dimensional vector of partial derivatives of /(A,r, r) by 

F = (d T1 f, ..., d Td f, d Xl f,..., d x J) T . (5) 


/(A,r,r) = 


/ 3^1 a—f t 2 d ) exp fy a itf + y 

is d -!(r) r Vi=l i=l 


Elements of F are called “standard monomials” in HGM. By Theorem 3 of [9j we have 


<9 r F = P r F, 


( 6 ) 


where the 2d x 2d matrix P r = (p tJ ), called the Pfaffian matrix, is of the form 



/ 2r 2 Ai + 1 

o 

P 

Ti A 

1 

O 

2 r 2 X d + 1 r d 

Td 

r 

r 2 Ti 

O 

2r 2 Ai+ 2 

1 


V O 

r 2 r d 

1 

2r 2 A d + 2 ) 


(7) 


with O denoting an off-diagonal block of 0’s and 1 denoting an off-diagonal block of l’s. The 
elements p t] of P r are expressed as 


d 

rpij = (2A;r 2 + 1)5^ + y pS j{k+d) (1 < i < d), 

k =i 

fP{i+d)j PSij T (2Ajr" T 2')8j/is r( p T ^ ^ $j(k+d) (1 F: i /X d), 

k^i 


for 1 < j < 2d , where S tJ denotes Kronecker’s delta. Given initial values for the elements of F at 
r = r 0 , we can apply a standard ODE solver to flSD for numerically evaluating F. 

For the initial values at a small r = ro > 0, we can use the following series expansion of the 
Fisher-Bingham integral my 


/(A, r,r) 


= r 


d -1 


Sd-1 


X ^ ' r ^\ a +p\ 

a,/3eNg 


(d- 2)!!nt 1 (2w + 2A-l)!! 2/? 

(d — 2 + 21 a | + 21/31)! !cv! (2/5)! 


( 8 ) 
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where No = {0,1,2 ,...} and for a multi-index a £ Ng we define 


d d d 

a\ = a!! = and |a| = otj. 

2—1 2—1 2 — 1 

By term by term differentiation of this series we can evaluate derivatives of /(A,r, r). For com¬ 
puting the initial values, we apply the following approximation: 

S' d _ir d+1 r i + 0(r d+3 ) (i = 1,..., d), (9) 

S d -ir d+1 + 0(r d+3 ) (i = l ,...,d). (10) 

By this approximation, we reduce the computational time for the initial values. However the 
accuracy of the result does not decrease at all as we will show in Section [H 

As r —* oo, the absolute values of /(A, r, r) and its derivatives become exponentially small, as 
we analyze the behavior in the next section. Hence we also consider the following vector 

Q = exp(-r 2 Ai - 7 ~|ti|) d T2 f ,..., d T J, ^d x J, d x J ,..., d x J^j . (11) 

Then from (jUJ) it is easy to obtain <9 r Q as 

d r Q = (D~ 1 d r D - (2rA, + |n|)/ 2d + DP r D- 1 ) Q, (12) 

where I 2 d is the identity matrix with size 2d and 

D = diag 1,4 1 ,-A 

\r r z ) 

The equation (TTTh is a refinement of the equation (21) in [9]. By Proposition 13. ll in the next section, 
each element of Q converges to some non-zero value when r goes to the infinity. This prevents the 
adaptive Runge-Kutta method from slowing down. 


df_ 

dri 

df_ 

d\i 


3 Laplace approximation close to the infinity 

In our implementation of HGM, we start from a small r = > 0 and numerically integrate F in 

([5]) up to r — 1 and then integrate Q in (ITT} toward r = oo. In order to assess the accuracy of 
Q for large r, we derive the asymptotic values of the elements of Q by the Laplace method. The 
Laplace approximation, including higher order terms, for the Fisher-Bingham integral itself was 
given in ra- However here we also need approximations for its derivatives, which were not given 
in [12J. Hence we give the approximations of the main terms of the Fisher-Bingham integral and 
its derivatives and a sketch of their proofs. 

We first consider the case of single largest Ai- We state the following result. 
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Proposition 3.1. Suppose 0 > Ai > A 2 > • • • > A d- Then, as r —> 00 


/(A,r,r) = 


n (d- 1)/2 


Ili 2 (Ai - Ai ) 1 / 2 
d\J{\,T,r) = r 2 /(Ai,ri,r)(l + o(l)), 


(e rT1 + e.- rTl ) exp r 2 Ai - 


T? 


£ ^ ~ ^ 


(1 + 0 ( 1 )). 


d\J(\,T,r) = 


Ti 


+ 


2(Aj-Ai)y 2(Ai - \j) 


/(A, t, r)(l + o(l)), (j — 2,... ,d) 


0 rr\ _ p—rri 


drj( A, r, r ) = A, r, r)(l + o(l)), 


e rn _p g-rn 


d Tj f(\r,r) = 


Ti 


2(A, - Aj) 


/(A,r,r)(l + o(l)), (j = 2,...,d). 


(13) 

(14) 

(15) 

(16) 
(17) 


Note that for Ti > 0, in (TT3l) e~ r ' Tl is exponentially smaller than e rT1 and it can be omitted. 
However we leave e~ rT1 there for consistency with the case of Ti = 0. Also we found that leaving 
e~ rTl in (fTSD greatly improves the approximation. 

We now give a rough proof of Proposition 13.11 . In the proof, the main contributions from the 
neighborhoods of maximal points are carefully evaluated, but the contributions from outside the 
neighborhoods are not bounded rigorously. Replacing ti by and integrating over S' a!_1 (l) can 
write 


/(A, t, r) = r d ~ 
d\J(\ ,r,r) = r d+1 
dr,/(A,T,r) = r d 


r / d d \ 

1 exp ( r 2 A itf + t V] Titi ) dt, 

Js*-HD \ iTt tT / 

/ t 2 exp ( r 2 V A itf IrV t+ ) dt, 

r ( d d \ 

/ t 3 exp ( r 2 Y" A+ 2 IrV ++ I dt. 

ds^-hi) V i =1 J 


(18) 

(19) 

( 20 ) 


For very large r 

t 2 (AR 2 + A 2 t 2 +-f A d t 2 d ), 1 = t\ +-1- 4 (21) 

takes its maximum value at two points t± = ±l,i 2 = • • • = td = 0. The main contributions to 
m-m come from neighborhoods of these two points (±1, 0,..., 0). The contribution from the 
complement of these two neighborhoods should be exponentially small as r —> 00 , although we do 
not give a detailed argument. We also have to consider the effect of r Yli=i TU- But it is of the 
order 0(r), whereas (1211) is of the order 0(r 2 ). Hence r ^2'!= 1 A+ only perturbs the maximizing 
values (±1, 0,..., 0) by the term of the order 0(l/r). Based on these considerations write 


t 


2 

1 


1 -t 


2 

2 


4 h = ±\Ji - ti - t 2 d = ±(! - ^(4 + • • • + 4 ), 


where |t 2 |,..., \td\ are small. As shown below, \ti\, i — 2,..., d, are of the order 0{l/r). We now 
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consider the neighborhood of (1, 0,..., 0). By completing the squares we have 


2—1 


A ^ 2 +r Yl Titi 
2—1 

d 

Ai + rTi + r" ((A i — Ai — — )t 2 + —ti) + o(l) 

1=2 
d 

= r 2 Ai + m + X! [( A * - Ai - + 


= r 


Ti 


T? 


i=2 

d 


= r 2 Ai + rri + ^ [(Ai - Ai)(rfj + 


2(Aj — Ax — ^) 4(Ai-Ai-§) 


] + o(l)- 


]+o(l) (22) 


i=2 


2(Ai — Ai) 4(A* — Ai) 


Furthermore around (1, 0,, 0) the volume element dt of the unit sphere S d ~ 1 ( 1) is approximately 
equal to the Lebesgue measure df 2 • • • dtd, with the error of the order t 2 + • • • + t 2 d . Hence by the 
change of variables 

Ui = rt u i = 2,...,d, 

the contribution to /(A, r, r) from the neighborhood of (1,0,..., 0) is evaluated as 

_2 


exp(r 2 Ai + TTi) / exp((Ai - Ai)(tq + 


Ti 


;) 2 - 


2(Aj — Ai) 4(Aj — Ai) 


)du\... dud 


= exp [ r 2 Ai + 


t. 


7r 


(d- 1)/2 


(23) 


^4(A,-A 1 )y nt 2 (Ai-Ai)V2' 

Similarly by changing the sign of T\ we can evaluate the contribution from the neighborhood 
of (—1, 0 ,..., 0) as 

d 


t; 


n (d- 1)/2 


exp|r2Al _ rTl _g 4(A __ A ^ nt ^_ Ai)i/2 . 

Adding (123|) and (124)1 we obtain (115)) . 

For d\ 1 f(X,r,r) and d Tl f(X,r,r), we can just put t\ = ±1 in (1T9)) and 
tions from two neighborhoods we obtain (114)) and (TT6)) . 

For d Xi f( A, r, r) and d T J( A, r, r), j > 2, we write 


(24) 

Adding contribu- 


u i 


tj — — — - [ U-; + 


Ti 


t 2 = — 

3 r 2 


2(A j — Ai) 2(Aj — Ai)/ 

2 


Ui + 


Ti 


Ti 


J 2(Aj-Ai) 2(Xj - Ai) 

and take the expectation with respect to a normal density. Then we obtain (115)) and (ITT)) . Although 
we did not give a detailed analysis of the remainder terms, we can show that the relative errors in 
(TT4D - (fT7)) are of the order 0(l/r). This completes the proof of Proposition 13.11 

A generalization of Proposition 13.II to the case that Ai = • • ■ = X m > A m+i > • • • > A^ is given 
in Appendix. We note that numerically HGM works fine even if some of the A’s are close to one 
another, because the Pfaffian system does not have a singular locus except at r = 0 and the main 
exponential order is the same in Proposition 13.11 and in Proposition IA.1I However when we want 
to check whether the ratio of HGM to the asymptotic value is close to one, then we have difficulty 
when some of the A’s are close to one another. 
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4 Numerical experiments 

In this section we describe our numerical experiments on the performance of HGM. The programs 
and the raw data of our numerical experiments are obtained at 

http://github.com/tkoyama-raaylO/ball-probability/ 

Our programs utilize the Gnu Scientific Library[6]. 

In our experiments we compute the initial values of r _ ^ +1 ^F at r = r 0 = 1.0 x 10 -6 by (J^J) and 
(Hop . The reason for multiplying F by r~^ d+1 ' is that the values of elements of F are too small at 
r 0 for floating point numbers. Then up to r = 1, we solve the differential equation ([6]) numerically. 
In our implementation, we utilize explicit embedded Runge-Kutta Prince-Dormand (8, 9) method 
and we set the accuracy to 1.0 x 10 -6 . In order to prevent the elements of F becoming too large, 
we re-scale the elements of F several times. Then at r = 1 we switch to Q in (HU) and solve (TT2]h 

Note that we can not take r = 0 as an initial point. The point r = 0 is in the singular locus 
of the differential equation (j6j) since the denominator of P r becomes zero at the point. Hence, 
numerical differential equation solvers can not compute the differentiation of F at the point. 

Our implementation computes the initial value of F by the approximations ([9]) and (ITU , which 
use only the first term of the series expansions. Hence, we have to take very small value for r in 
order to reduce the error. 

Each component of vector F takes very small value at r o = ICE 6 . We deal with this problem 
by storing not the value of F itself but the product of F and a large constant in double precision 
type array. Related to this problem, there is another problem that the value of each component 
of F increases rapidly when we solve ordinary differential equation (J6]) numerically. We multiply 
vector F by a small constant when a component of F becomes larger than a fixed value. By this 
way, our implementation prevents values in double precision type array becoming too large. 

Our first experiment is for d = 3 and the following parameter values 


i.e., 


(7 1 = 3.00, 

/i, = 1.00, 


Ai = -0.0555556, 

Ti = 0.111111, 


(72 = 2.00, 
p 2 = 0.50, 


A 2 = -0.125, 
t 2 = 0.125, 


cr 3 = 1.00, 

p 3 = 0.25, (25) 


A 3 = -0.5, 
r 3 = 0.25. 


By HGM we compute G(r). We show its graph in Figured] to confirm that our implementation 
correctly calculated the asymptotic behavior as G(r) —> 1 as r —> oo. 

For this example, we also check the accuracy by computing the ratios of /(A,r, r) and the 
elements of F to their asymptotic expressions in Proposition 13.11 The left figure of Figure [2] shows 
the ratio of /(A, r, r) to its asymptotic expression and the right figure shows the ratios of elements 
of F to their asymptotic expressions. Note that the value of the ratio corresponding to df/dXi 
is very close to that of df jdT{ so that the triangles overlap with the circles. We see that the 
numerical integration involved in HGM, starting from a small r 0 , is remarkably accurate, so that 
the ratios numerically converge to 1 as r —>■ oo. 

In our second example we consider diagonal matrices FA 1 ) and FA 2 ) with diagonal elements 

E) W = ^ < d), (26) 
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(Fisher-Bingham lntegral)/(Laprace Approximation) 



Figure 2: Ratios to the Laplace approximations 


and 



2(d + 2)(d + 3) 
k(k + 1 )(k + 2){k + 3) 


(1 < k < d), 


respectively. These weights are considered for cumulative chi-square statistics in [7], Let 


h (1) = o, 

/i (2) = (0 0.01 0.02 


0.01 x (d- 1 )) T . 


( 27 ) 










For each dimension d, we computed the probability P (10 6 < ||X|| < 40.0) and measured the 
computational times in seconds. We considered the following four patterns of parameters: 

(£ (1) ,/i (1) ), (£ (1) ,/i (2) ), 

(£ (2) ,/i (1) ), (£ (2) ,/i (2) ). 

The experimental results are shown in Tablc|Tl 1 —p stands for the values 1 —P (10~ 6 < ||X|| < 40.0) 
are generally accurate to 10~ 8 . 


Table 1: Accuracy and computational times for E^ 1 ) anc [ y;( 2 ) 


dimension 

h 

1 — p 

= 0 0 
times(s) 1 — p times(s) 

h 

1 — p 

= 0 

times(s 

EW 

) 1 — p times(s) 

10 

1.60e-08 

0.03 

1.60e-08 

0.03 

1.60e-08 

0.11 

2.10e-09 

0.11 

11 

1.76e-08 

0.03 

1.57e-08 

0.04 

1.76e-08 

0.12 

1.56e-09 

0.14 

12 

1.61e-08 

0.04 

1.15e-08 

0.04 

1.61e-08 

0.16 

9.59e-10 

0.17 

13 

1.81e-08 

0.04 

1.05e-08 

0.04 

1.80e-08 

0.20 

7.90e-10 

0.19 

14 

2.02e-08 

0.04 

9.95e-09 

0.05 

2.02e-08 

0.24 

6.94e-10 

0.25 

15 

2.34e-08 

0.04 

9.58e-09 

0.06 

2.34e-08 

0.30 

6.44e-10 

0.30 

16 

2.77e-08 

0.06 

9.73e-09 

0.07 

2.77e-08 

0.36 

2.89e-10 

0.36 

17 

3.40e-08 

0.07 

4.85e-09 

0.08 

3.40e-08 

0.41 

2.74e-10 

0.42 

18 

1.89e-08 

0.08 

4.62e-09 

0.08 

1.89e-08 

0.49 

2.82e-10 

0.52 

19 

2.08e-08 

0.08 

4.40e-09 

0.10 

2.09e-08 

0.56 

4.05e-10 

0.57 

20 

2.33e-08 

0.10 

4.32e-09 

0.11 

2.41e-08 

0.65 

1.13e-09 

0.65 


As the radius r increases or the dimension d of the sphere increases, our implementation takes 
long time to evaluate. Table 1 shows that the computational complexity also depends on the values 
of A. However we do not know what value of A makes the computational time worse. 

As our third example we consider how our method works for large dimension. Corresponding to 
the asymptotic null distribution of Anderson-Darling statistic, which is an infinite sum of weighted 
X 2 variables, consider the weights 

al = WTTy , “ = 0 {1 - k - 


Here we truncate the infinite series at d. We computed the probability and measured its com¬ 
putational time. We fixed the radius as r = 20.0. The results on the computational time are 
shown in Table [2] and its figure. Even for d = 100, our method is accurate and fast enough to 
be practical. This is a remarkable progress since the implementation of HGM in [9] can compute 
only up to dimension d = 8. The key idea for this progress are the simple approximation of the 
initial values (J9]) and (TTOi) for HGM and the refined differential equation (fl2|) based on the Laplace 
approximation. 

The computational bottleneck of HGM is the computation of P r F in each step of solving the 
ODE. By the form of the matrix P r , the number of additions in each step increases in order 0(d 2 ). 
We guess this is a reason that growth of computational times in the figure of Table [2] seems to be 
in the order 0(d 2 ). 
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Table 2: Computational times for Anderson-Darling statistic 
dim 1 — p time(s) 


30 

5.70e-08 

1.03 




35 

3.76e-08 

1.59 




40 

4.85e-08 

2.36 


OD " 


45 

6.13e-08 

3.30 



O 

50 

8.97e-08 

4.42 

CD 

E 

LO 

CM 

O 

55 

5.29e-08 

5.94 

IB 

c 

O 

CM 

o 

60 

7.91e-08 

7.56 


LO 

o 

65 

6.28e-08 

9.69 

Q. 

E 


o 

70 

1.02e-07 

12.05 

O 


o 

o 

75 

6.77e-08 

14.63 


LO - 

o 

n ° 

80 

7.22e-08 

17.81 


o - 

> 

o 

o 

o 

o 

o 

o 

o 

85 

6.25e-08 

21.33 



1 1 1 1 1 

20 40 60 80 100 

90 

5.64e-08 

25.10 




95 

5.21e-08 

29.54 




100 

4.90e-08 

35.05 


Graph of computational times 


As our fourth example we consider the case where S is the identity matrix and /i = 0. In this 
case, the Fisher-Bingham integral can be written by the density function of ^-distribution, and 
we have 


f(r) 


2n d / 2 

nd/2) 


r d-l e -r 2 /2 


Table \3\ shows the result for /(1.0) by HGM and difference 1 e r2//2 — /(1.0) for each 

dimension. 


Table 3: Comparison to ^-distribution: T, — I and /r = 0. 


dim 

hgm 

exact—hgm 

3 

7.621888 

1.35e-06 

4 

11.972435 

7.09e-07 

5 

15.963247 

4.85e-07 

6 

18.806257 

3.70e-07 

7 

20.060008 

3.34e-07 

8 

19.693866 

3.13e-07 

9 

18.005821 

2.88e-07 

10 

15.467527 

2.47e-07 


As our fifth example we consider the case where 
^ _ / i 1 1 1 1 1 

In this case, the ball probability (JT|) equals to 


\/2 ~n y/2n) 


), 


H = 0 (d — 2 n). 


P [ -X 2 + - X 2 + -X 2 + -Xl + -Xl + -X 2 + 
' 2 1 2 2 4 3 4 4 6 5 6 6 


+ + < r2 
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where Xi ,..., X 2n are independent and identically distributed with the standard normal distri¬ 
bution. Since the distribution of A (X| fc _ 1 + A| fc ) is the exponential distribution with the rate 
parameter k, the above probability is equal to (1 — e~ r2 ) n [5j p.21]. The second column in Table [4] 
shows the result of HGM for the ball probability at r = 1.0. The third column shows the difference 
between HGM and the exact value. 

Table 4: Comparison at specific parameters. 


dim 

hgm 

exact—hgm 

6 

0.252580 

4.97e-09 

8 

0.159661 

2.54e-09 

10 

0.100925 

1.61e-09 

12 

0.063797 

1.03e-09 

14 

0.040327 

8.16e-10 

16 

0.025492 

7.07e-10 

18 

0.016114 

3.04e-10 

20 

0.010186 

2.37e-10 


5 Summary and discussion 

In this paper we applied HGM for computing distribution function of a weighted sum of inde¬ 
pendent noncentral chi-square random variables. We found that our method is numerically both 
accurate and fast, after we implemented the following ideas. First, during the application of Runge- 
Kutta method, we re-scaled the vector F in ([5]) as needed to keep its elements within the precision 
for floating point numbers. Also we divided the interval for integration into (0,1] and [1, oo) and 
switched from F to Q in (1111) in view of the asymptotic values for Q. Our experience in this 
paper shows that re-scaling of the standard monomials is important in numerical implementation 
of HGM. 

In our implementation, the numerical integration starts from a small r = ro > 0 and the 
integration proceeds to r = oo. On the other hand, we have asymptotic results for large r in Section 
[3] Then we might consider reversing the direction of integration and start with initial values at 
very large r. We may call the former the “forward integration” and the latter the “backward 
integration”. However we found that the backward integration is not numerically stable. Hence 
the asymptotic values can not be used as initial values. In this paper we used the asymptotic 
values just for checking the accuracy HGM in the forward direction. 

It is an interesting question, whether the asymptotic values can be used to adjust the values 
of the forward integration. We may look at the difference between F by forward HGM and its 
asymptotic value for very large r and use the difference to adjust F at intermediate values of r. 
However it is not clear how this adjustment can be implemented. 

A A general form of Proposition 13.1 

In Proposition 13.II we assumed Ai > A 2 . In this appendix we state the following proposition for the 
general case Ai = • • • = A m > A m+ i without a proof. For this case, the integrand for the Fisher- 
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Bingham integral takes its maximum on the (m — l)-dimensional sphere S m ~ 1 ( 1), rather than on 
a hnite number of points. However by appropriate choice of coordinates and by multiplication of 
the volume Vol(S' m_ 1 (l)), the derivation of Proposition IA.1I is basically the same as Proposition 

o 

Proposition A.l. Assume that 


0 > Ai — • — A m > \m+i + ■ ■ ■ + Ad. 

If 0 = T\ = ■ ■ ■ = T m , then as r —$■ oo, 


/(A,r,r) = r m 1 S m -i exp I r 2 Ai - ^ 


?r ( d-m)/2 


i=m -\-1 


nt m (Ai-A 4 )V2 


( 1 + 0 ( 1 )), 


d\J(\,r,r) = —/(A,r,r)(l + o(l)), j < m, 


m 


d T J(\r,r) = 0 , j<m, 


d T J(\r,r) = 


To 


2(A j — Ai) 


/(A,r,r)(l + o(l)), j>m, 


d\J(\,r,r) = 


+ 


r. 


/(A,r,r)(l + o(l)), j>m. 


,2(Ai — Xj) 4(Aj — A 1 ) 2 / 

If (ti, ..., T m ) ± (0,..., 0), de/me 7 = (rf H-b r,^) 1/2 . Then, as r ^ 00 , 


/(A,r,r)=exp r 2 Ai+ry- ^ 


z=m+l 


4(Aj - Ai) / V 7 


2 r 


( m — 1)/2 


7T 


(d-l)/2 


n, =m (Al - Ai) 1 / 2 


(l + o(l)), 


<9 Ti /(A,r,r) =r-4/(A,r,r)(l + o(l)), t.,- 0 0, j < m, 

7 + 

fty/(A, r, r) = r 2 -+/(A, r, r)(l + o(l)), t,- ^ 0, j < m, 

d Tj f{ A, r, r) = 0 , Tj = 0 , j < m, 

^Aj/(A,t, r) = —/(A, r, r)(l + o(l)), r,- = 0, j < m, 


7 


^/(A,r,r) = 
d\J(X,T,r) = 


2 (Aj - A0 

1 


/(A,r,r)(l + o(l)), j > m, 


+ 


2(A 1 -A i ) 4(Aj — Ai) 2 


/(A, t, r)(l + o(l)), j>m. 
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